Mathematical Biosciences and Engineering

In many natural dynamical systems, tipping points emerge where slowly changing external factors induce shifts from one state to another, known as critical transitions or bifurcations. In epidemic systems, detecting these transitions is important for outbreak preparedness. Various statistical indicators like lag-1 autocorrelation and variance have been employed as early warning signals (EWS) of critical transitions by showing how the dynamics of the system exhibit critical slowing down (CSD) near the transition. More recently, deep learning-based algorithms have been developed that attempt to more accurately detect and classify bifurcations; however, they are limited by the systems simulated in the training dataset, the length of the time series, and have not been rigorously evaluated on real-world epidemic data spanning multiple geographic regions. In this paper, we evaluate existing deep learning-based bifurcation prediction models, on a CDC-curated dataset containing weekly influenza hospital admissions across all fifty U.S. states and the nation as a whole. We develop a pipeline for data preparation, reproduction number (Rt) estimation using the EpyEstim package (a Python implementation of the EpiEstim framework [4]), and bifurcation point identification. We systematically evaluate the models under several windowing strategies, including rolling, expanding, and fixed-length windows, and report accuracy, precision, recall, F1 score, and specificity. Our best-performing configuration (expanding window, length 100) achieves an accuracy of 90.09%, precision of 72.86%, recall of 96.23%, and F1 score of 82.93%. We additionally discuss Markov Regime Switching models as an alternative framework for detecting epidemic transitions. Our findings indicate that these models can generalize to real-world, state-level influenza data and are capable of online/real-time prediction.

We develop an age- and immunity-structured COVID-19 model with setting-specific contacts, calibrated to Ontario incidence and hospital admissions before Omicron. We simulate counterfactual vaccination strategies and contact restrictions in schools, workplaces, and other community settings, evaluating projected severe-case demand against an illustrative healthcare-capacity benchmark. Vaccination combined with workplace and community contact restrictions keeps simulated demand within capacity. In contrast, relaxing restrictions in either setting produces large surges and repeated capacity threats, even when school restrictions remain. School-only changes produce smaller effects on severe disease. These findings highlight the value of vaccination and setting-specific contact restrictions that prioritize workplace and community settings for pathogens with adult-concentrated severe outcomes.

Mechanistic models of infectious disease transmission provide a framework for understanding and predicting disease dynamics, but statistical inference is often computationally prohibitive. This challenge is especially pronounced for models with high-dimensional structures and many unknown parameters, where trajectory matching may require thousands of evaluations of an expensive likelihood function. We developed a surrogate-assisted optimization framework that uses Bayesian optimization to reduce the number of evaluations of the underlying mechanistic model required for parameter inference. We evaluated the approach using simulated fitting problems based on age-structured susceptible-exposed-infected-recovered models with 4, 8, and 16 age classes and a 210-compartment age-structured rotavirus transmission model. Applying surrogate-assisted optimization substantially reduced the number of likelihood evaluations required to identify good parameter estimates, with larger improvements as model dimensionality increased. In the rotavirus model, the surrogate Bayesian optimization approach reached the benchmark after approximately 0.7 hours, compared with approximately 9.4 hours using traditional Nelder-Mead optimization. These results demonstrate that surrogate-assisted optimization can substantially reduce the computational burden of trajectory matching for complex mechanistic transmission models, particularly for higher-dimensional inference problems.

Blooms of the toxic dinoflagellate Ostreopsis cf. ovata can cause respiratory irritation in beach users, yet their environmental drivers and temporal dynamics remain poorly understood. Here, we use a dynamical-systems approach to investigate how environmental conditions shape bloom dynamics from long-term time series of algal abundance and abiotic variables. Using the global modelling framework GPoM, we reconstruct nonlinear ordinary differential equation systems and identify candidate models capable of reproducing key features of the observed dynamics, including chaotic behaviour. Topological and statistical analyses reduce the candidate models to two systems that best capture the observed distributions and dynamical patterns. These models reveal distinct combinations of environmental conditions associated with bloom occurrence and amplitude, providing a framework for interpreting bloom trajectories and a potential basis for short-term prediction.

Automated surveillance systems can prove to be extremely useful in detecting outbreaks of respiratory diseases that need further epidemiological investigation. In order for the system to be effective, careful attention has to be paid to the baseline used, absence of the observations, and the use of proper statistical methods. As a solution, the paper describes EpiAlert, a weekly surveillance system based on CDC FluView ILINet. The system is able to transform the number of ILI visits and the total number of visits into a ILI prevalence for the region, create the first reference from the three previous years using a +/- 2-week interval, transform the rate value, and apply two different detectors - CUSUM and EWMA. CUSUM is set at k=0.5 and h=5, whereas EWMA uses {lambda}=0.2 and H=2.7 values. Among the 55 districts, the system gives 45,138 pieces of information about the disease by transforming 44,297 pieces of data into alerts. It is found that there are 2,820 alerts in total (out of which 2,711 alerts used CUSUM, 402 used EWMA, and 293 used both detectors). Additionally, the experiments show that the sensitivity of the device depends greatly on the detector. The latest statistics can be tracked in the report for week 37 of 2026 that recently concluded in West.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

Dynamic systems often evolve faster than reliable measurements can be collected. More frequent indirect observations may fill these gaps, but their noise and differing scales complicate estimation of the underlying state. We developed Bayesian Anchored Latent Learning (BALL) to combine sparse reference measurements with denser records on an interpretable scale. A bidirectional transformer teacher reconstructs trajectories from complete training histories, and a forward-only transformer student learns from those reconstructions and observed measurements. Across 30 simulations, BALL reduced composite trajectory error by 7% relative to the same transformer trained directly from measurements, achieving lower error across all datasets. In an applied example using psychiatric records from 1,609 patients, BALL improved the prediction of withheld questionnaire scores compared with direct transformer training at 21- and 28-day measurement schedules. Requests guided by estimated uncertainty further improved prediction. BALL connects retrospective learning with prospective state estimation and provides a framework for directing additional measurement.

Fractional differential equations are increasingly used to represent temporal memory in tumour immune and pharmacokinetic models, but changing the derivative order without a consistent time scale can make comparisons difficult to interpret. We study a dimensionless, commensurate Caputo system for tumour, immune, normal cell, and drug states. We establish positivity, boundedness, the tumour invasion threshold, and local stability conditions. We then use a fractional Adams Bashforth Moulton predictor corrector method with an explicit mesh convergence study. Two conclusions correct potentially misleading interpretations of the model. First, the equilibrium locations are independent of the fractional order: decreasing the order changes transient time scales but does not by itself alter the model's asymptotic states. Second, after a unit bolus the drug equation has the exact solution C({tau})=E(-K{tau}). Its maximum is therefore the initial concentration, not a delayed peak. For 0<<1, the algebraic tail makes the infinite-horizon area under the concentration curve diverge. In the illustrative coupled simulations, the apparent late tumour reduction at small is consequently driven by persistent modelled drug exposure and must not be read as evidence of superior therapy. All parameters are dimensionless and illustrative; no patient data are fitted. Source code and generated tables are included to make every numerical statement reproducible.

Pseudoprogression in response to immune checkpoint inhibitor (ICI) therapy can be a complex challenge during treatment response evaluation. Erroneous classification of progression may result in stopping effective treatment for patients, while delayed confirmation of progression can prolong an ineffective treatment course. We develop a mechanistic mathematical model of tumor-immune dynamics under ICI treatment in which observed tumor volume includes contributions from both tumor burden and immune infiltration. We show analytically how pseudoprogression can arise from increased immune infiltration and proliferation as well as continued tumor growth during the ``ramp-up" phase. The model also incorporates circulating tumor DNA (ctDNA) dynamics, which we demonstrate are able to closely fit existing longitudinal tumor and ctDNA data and extract realistic parameter ranges to generate representative in silico data. Using model-generated cohorts, we train a random forest (RF) classifier using longitudinal tumor and ctDNA measurements to distinguish pseudoprogression and true progression, and demonstrate that it improves response classification relative to the Response Evaluation Criteria in Solid Tumors (RECIST 1.1) and provides classification earlier than the immunotherapy-specific iRECIST, even in the presence of tumor measurement noise. This work predicts that paired longitudinal tumor-ctDNA measurements may enable earlier discrimination of true progression from pseudoprogression, and motivates collection of imaging and ctDNA data synchronously during ICI therapy, especially in the setting of suspected disease progression.

The 2025 measles outbreak in South Carolina, the largest in the United States in over two decades, produced 997 confirmed cases, predominantly among school-aged children in an under-vaccinated community in Spartanburg County. To support the South Carolina Department of Public Health (DPH) measles response, we developed a stochastic agent-based model to forecast outbreak trajectory and evaluate intervention strategies in real time, integrating transmission across three spatial scales and grounded in geocoded households, school enrollment records, and census-derived demographic data. The model predicted a median of 1,020 total infections (90% prediction interval: 840-1,273). Counterfactual simulations without active contact tracing and quarantine projected a median of 2,432 infections (90% prediction interval: 1,867-3,185), approximately 138% more than the model-predicted baseline (1,020 cases). Model fidelity required a highly localized between-school contact structure; broader contact assumptions substantially overpredicted both case counts and geographic spread. Simulated regional introductions confirmed that outbreaks arose only near connected clusters of undervaccination, while well-vaccinated areas remained contained. Model versions deployed prospectively at weeks 9, 17, and 34 of the response informed resource allocation decisions and provided operational evidence for the effectiveness of contact tracing and quarantine strategies. Critically, the model reproduced geographic spread, number of affected schools, and breakthrough infection counts not used in fitting, validating its structural fidelity and establishing a replicable framework for jurisdictions maintaining school-enrollment and vaccination records.

Crickets are gaining attention as sustainable protein sources because of their high nutritional value and low environmental impact. Growth models are essential tools for optimising production by predicting growth characteristics. However, previous studies have not compared multiple non-linear growth models for crickets or explicitly considered environmental factors. We analysed the growth patterns of the two-spotted cricket (Gryllus bimaculatus) and the Asian cricket (Teleogryllus occipitalis) at multiple temperatures (25.0, 27.5, 30.0, 32.5 and 35.0 {degrees}C) using three non-linear models (Gompertz model, logistic model, and von Bertalanffy model) and combined the predicted body mass with observed survival rates to predict the optimal harvest time and temperature. All three models fitted the data well (R2 [&ge;] 0.97 in both species). The asymptotic body mass was highest at 30.0 {degrees}C, and the growth rate constant increased with temperature, reaching its maximum at 35.0 {degrees}C except for the logistic model in G. bimaculatus (32.5 {degrees}C). In addition, the predicted value of harvest, calculated by combining predicted body mass and observed survival probability, was highest at 30.0 {degrees}C in both species, reaching 13.82 kg/m3 at 38 days in G. bimaculatus and 7.12 kg/m3 at 62 days in T. occipitalis. These results suggest that 30.0 {degrees}C is the most suitable constant-temperature condition for maximising predicted productivity, although the optimal harvest timing differs between species.

The emergence of Highly Pathogenic Avian Influenza (HPAI 2.3.4.4b) in dairy herds in 2024 across 19 states in the United States of America has raised concerns regarding the potential national and global zoonotic impact. All recent modeling efforts implemented homogenous cattle-to-cattle (both intra and inter-herd) transmission models which did not capture the real-world heterogeneity in mixing of animals, individual variations in susceptibility and infectiousness and clinical incidences across pens and lactation groups. The aim of this study was to develop a heterogenous transmission model to estimate the epidemiological parameters for intra-herd HPAI transmission on Californian dairies. We developed a validated stochastic agent-based model to estimate the epidemiological parameters for intra-herd HPAI transmission on California dairies. The hierarchical agent-based model also parameterized stochastic cattle movements within-herd to simulate real dairy management practices. A cow-level SEIR transmission approach was assumed during the outbreak. A novel Bayesian Optimizer with Gaussian Process (BO-GP) was fitted to the agent-based model for validation which converged within 25-40 iterations (out of 100 per farm) with minimal loss over two distinct error metrics, namely, Poisson loss function and temporal distance metric. Our optimized simulations estimated an average R0 was around 10.7-10.8 across all farms within the first 15 days of observed outbreak on four dairy farms with a mean effective transmission rate of 4.5% per contact between susceptible and infectious cows within each pen. Our model demonstrated that the movement of cows between pens ensured localized clusters of outbreaks within the sub-herds (pen population) that prolonged the overall outbreak within farms. We estimated the total duration of infection between 14.5 and 28 days, which is higher than the estimates from the homogenous models. With an integrated hierarchical agent-based model combined with Bayesian approximation, we produced actionable insights on the epidemiology of intra-farm spread of HPAI within cow herds, thereby guiding both future model development and applied disease control strategy.

We propose a stochastic epidemic model to describe the spread of avian influenza in a network of poultry farms. We consider three main paths of infection transmission between poultry farms: from nearby farms, mainly through airborne diffusion; from farms belonging to the same company, through shared veterinarians, forage providers and similar; from undetected small farms or wildlife. From the parameter estimation, based on a modified Expectation-Maximization algorithm, we infer that approximately 63% of the farms were infected from nearby ones, with an infection force declining with distance; of these, more than one third belonged to the same company of the estimated infector. About 20% were infected from premises belonging to the same company but farther away than the distance threshold of 2 km; the remaining ones from wildlife or unidentified sources. These estimates have been validated by a comparison with genetic data, available for a subset of the farms: the genetic distance between two farms identified, with high probability, as an infector-infectee pair is much lower than between random pairs. These results may help in the implementation of tailored prevention and control measures for future outbreaks.

In this study, we develop a mathematical model of the hypothalamic pituitary adrenal (HPA) axis and systematically investigate the model complexity required to reproduce observed hormone dynamics. The model is fitted to ACTH and cortisol measurements from five healthy volunteers reported in the literature. We show that simultaneous reproduction of both hormone profiles requires a delay differential equation model with at least five state variables among the model architectures considered. Following fitting the volunteer data, we investigate the response of the model to sustained stress using bifurcation analysis. The coupled system exhibits bifurcations over a range of stress intensities. This analysis reveals bifurcations whose critical stress thresholds depend strongly on the modulatory state of the system. These results provide a mathematical framework for investigating transitions between physiological and dysregulated HPA axis dynamics, which may be related to states of health and disease.

Oncogene amplification on circular extrachromosomal DNA (ecDNA) has been linked to poor prognosis and higher treatment resistance in multiple types of human cancer. ecDNA are mobile genetic elements lacking centromeres that are partitioned unevenly into daughter cells at mitosis. While random segregation of ecDNA contributes to gene-copy-number heterogeneity among tumour cells, how ecDNA contribute to phenotypic heterogeneity, via mRNA and proteins, and how this affects targeted therapy outcomes is still not understood. In fact, cancer cell populations with ecDNA show remarkable heterogeneity in mRNA concentrations, which is not explained by gene-copy-number heterogeneity alone. In this paper, we use a stochastic model of gene transcription to mathematically model mRNA production and decay in cancer cells with ecDNA, extending a method that has been used to study phenotypic heterogeneity in eurokaryotic and prokaryotic cells. Using our model, we explore the space of possible mRNA distributions from a given ecDNA distribution, under different model parameter regimes. Our model is able to produce realistic mRNA distributions and we use it to estimate transcriptional parameters from single-cell measurements of ecDNA and mRNA counts in cell-line data. We discover a remarkable similarity in estimated parameters across cell lines and amplification types, indicating similar processes underlying transcription in ecDNA and chromosomal amplifications. We also model genotoxic therapy and use our model to predict the mRNA composition of persisting cells, showing the usefulness of such a modelling framework to understanding the origins of treatment resistance in ecDNA-containing cancers.

An increasing number of researchers, educators, and at-home growers are turning to hydroponics as an alternative method of food production, diversifying and strengthening the food supply chain. This has led to a substantial increase in the amount of academic work involving hydroponic production of plants for a variety of purposes, including but not limited to a more controlled and customizable environment than conventional production methods, utilizing hydroponics as a method of wastewater treatment, and investigating the vertical production of plants on a much smaller footprint than traditional commercial-scale production. Because so many unique stakeholders are interested in hydroponics for a variety of reasons, academic studies use a wide variety of hydroponic systems with varying levels of technology and automation, making comparisons across studies and systems challenging. A few institutions have sophisticated plant growth and research facilities, like those used by commercial companies, but these systems are often prohibitively expensive, costing upwards of $10,000 USD for a very limited amount of growing space for experimental trials. To enable more researchers, educators, and citizen scientists to conduct experiments and easily compare their results, this study details three pilot-scale hydroponic production systems (Kratky, NFT, and DWC), utilizing entirely off-the-shelf components, for under $5,000 USD each, that are capable of maintaining the desired consistent temperature and relative humidity levels, and have space for 216, 108 (6 irrigation loops), and 64 (16 irrigation bins) plants, respectively. The material lists and prices are provided along with assembly instructions and general guidelines for each hydroponic system. The testing and optimization process of environmental conditions within each growth system is also detailed, with suggested customizations for potential experimental setups. The authors developed these systems during their Masters and Doctoral degree programs, enabling improved hydroponic research at the University of Illinois Urbana-Champaign. They hope that other aspiring hydroponic researchers can benefit from this work in the development and operation of their hydroponic systems, increasing research efficiency and collaboration opportunities in this field of study.

Undersea and high-altitude operations regularly put operators at risk for decompression sickness (DCS), which is attributed to bubble formation and growth during decompression. DCS mitigation strategies are validated with empirical testing, often conducted in hyperbaric chambers, which provide precise control of the pressure-time profile. A degree of variability for study outcome is anticipated in both decompression sickness (DCS) and venous gas emboli (VGE) endpoints. Experimentation in a chamber offers precision, but is not always practical or operationally relevant, leading researchers to use open water for the dive exposure. Dive profile consistency may vary by diver skill level and environmental conditions, potentially introducing dive profile deviations, which may be a factor of individual variability. There does not exist a quantitative framework for evaluating dive profile consistency and quantifying deviations in open water studies. In this work, we evaluate existing dive profile metrics and propose a workflow to estimate dive profile deviations using the Buhlmann dissolved gas phase model. The dive profile metrics included depth and time-based calculations, as well as surface gradient factor, degree-of-conservatism and the root mean square error of the Buhlmann tissue compartments. We developed a set of simulated dive profiles and obtained a set of 40 open water dives previously collected by certified divers, both of which were used to evaluate performance of the dive profile metrics. Depth and time-based metrics offered sensitivity to varying depths and time durations of the dive profiles, respectively. For the field dives, we found the Buhlmann root mean square error metric performed best at identifying the deviation profile, followed by pressure root time. These evaluations demonstrate the utility and limitations of previously published dive profile metrics and emphasize the importance of accounting for dive profile consistency in open water decompression trials.

Multiview data analysis has become increasingly important for integrating heterogeneous information from multiple data sources. Cooperative learning provides a flexible framework for multiview integration by bridging early and late fusion strategies. However, its use of squared loss makes it sensitive to outliers. In this study, we propose robust cooperative learning based on the soft-IPOD framework. The proposed method introduces an observation-specific outlier absorption vector into the cooperative learning objective, thereby reducing the influence of outlying observations while preserving information sharing across views. The resulting optimization problem can be formulated as a weighted Lasso-type problem using an augmented design matrix and can be naturally extended to more than two data views. Numerical simulations and a real-data application to labor onset prediction demonstrate the predictive performance of the proposed method in the presence of outliers.

Nitrogen dioxide (NO2) is a widespread air pollutant associated with the risk and prevalence of various respiratory infectious diseases. Epidemiological studies have reported associations between NO2 exposure and COVID-19 incidence, and biological evidence indicates that NO2 exposure increases vulnerability to infection. Motivated by these findings and by the scarcity of mechanistic epidemic models that incorporate air pollutant data, here we introduce an NO2- and mobility-informed metapopulation model and study the first COVID-19 wave in Oklahoma. County-level NO2 concentrations enter the infection rates through a first-order parameterization, alongside high-resolution mobility patterns capturing heterogeneous inter-county connectivity. Using numerical simulations, we demonstrate that higher environmental sensitivity produces earlier and more synchronized outbreaks across counties, while lower sensitivity yields delayed, asynchronous trajectories. When fitted to county-level COVID-19 incidence data, our model outperformed the classic SIR model (fitted independently in each county) in over 75% of counties in Oklahoma across all standard error metrics. The improvements were even more pronounced in forecasting tasks across epidemic stages, with county-level win rates against the SIR model of 59.74% (surge), 87.01% (peak), 80.52% (post-peak), and 88.31% (decline) by nRMSE. To isolate the contribution of NO2, we compared our model against a mobility-only model with NO2 removed and homogeneous infection rates across counties. The NO2 term improved county-level forecasts during the surge (94.81% of counties by nRMSE), post-peak (88.31%), and peak (62.34%) stages, though not during the early and decline stages, suggesting that NO2 information can generate better forecasts than mobility alone during periods of high incidence. Our results demonstrate that environmental exposure can enter the transmission structure of a mechanistic epidemic model directly and that doing so yields measurable forecasting gains.

Interception of moving objects is a universal pattern of animal behavior with many industrial applications. Inspired by the behavior of the characid fish Brycon guatemalensis, which catches fruits falling onto the surface of a river, we study one-shot (i.e., open-loop) interceptions using both stochastic optimal control theory and reinforcement learning. In the first instance, an agent actuates a control (wait time before launch and heading angle) to meet its target, given a noisy estimate of the target initial state. In a minimal setting, where launch velocity and actuated control are constants, we show that the optimal wait time is attenuated by noise in both the target measurement and control actuation, as well as agent energy conservation. We then derive an upper bound for the heading angle variance which scales as the ratio of the non-vanishing optimal interception error to the initial planar separation. Complementing our control framework, our deep reinforcement learning approach shows that the neural network representing the fish first learns the value function for the interception, then learns how to aim, and finally learns the optimal wait time. Beyond the specific problem at hand, our results may also have applications to problems with little time for feedback or corrections such as perching, landing, and docking.